In [5]:
plotSets = ['TransvPlane', 'TransvPsAngles', 'TransvPsMomenta', 'LongPsT']
plotDefs = bd.set_plot_defs_from_distrs([targetDistr], setNames=plotSets)
plotDefs[6]['lims2'] = [0., 100.]
_ = bd.plot_distr([targetDistr], plotDefs)
import os
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import BeamDynamics as bd
import SimulationData as sd
import RFTrackTools as rfttools
import copy
import json
/home/tia/Repos/GIT_PSIPositronProduction/BdPyTools/BeamDynamics.py:11: UserWarning: ROOT module not available.
warnings.warn('ROOT module not available.')
from importlib import reload
reload(bd)
reload(rfttools)
<module 'RFTrackTools' from '/home/tia/Repos/GIT_PSIPositronProduction/RunningSimulations/RFTrack/RFTrackTools.py'>
# %matplotlib inline
# %matplotlib notebook
%matplotlib widget
plt.rcParams['figure.figsize'] = [9.6, 6.4]
defaultColorCycle = plt.rcParams["axes.prop_cycle"].by_key()['color']
# plotFont = {
# 'family' : 'sans-serif',
# 'weight' : 'normal',
# 'size' : 12
# }
# matplotlib.rc('font', **plotFont)
# plt.rc('legend', fontsize=10)
DISTR_REL_PATH = '../../Data/Geant4/Pcubed/P3dist_baseline_p.ini'
targetDistr, _ = bd.convert_astra_to_standard_df(DISTR_REL_PATH, delim_whitespace=False, filterSpecsSelector=None, zProjection=None)
targetDistr.describe()
| Ekin | Q | betaRel | gammaRel | pdgId | px | py | pz | t | x | xp | y | yp | z | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| count | 137614.000000 | 1.376140e+05 | 137614.000000 | 137614.000000 | 137614.0 | 137614.000000 | 137614.000000 | 137614.000000 | 137614.000000 | 137614.000000 | 137614.000000 | 137614.000000 | 137614.000000 | 137614.0 |
| mean | 49.268378 | 1.000000e-09 | 0.995805 | 97.415810 | -11.0 | 0.057750 | 0.002026 | 47.963867 | 0.092658 | -0.007513 | 2.115502 | 0.002042 | 0.335502 | 0.0 |
| std | 118.740549 | 6.245236e-23 | 0.016748 | 232.369458 | 0.0 | 7.121371 | 7.055557 | 119.063066 | 0.006691 | 1.112848 | 448.462702 | 1.113585 | 447.292326 | 0.0 |
| min | 0.101143 | 1.000000e-09 | 0.550595 | 1.197931 | -11.0 | -87.036454 | -131.671848 | 0.000602 | 0.000000 | -13.017345 | -1570.004841 | -15.773611 | -1570.700587 | 0.0 |
| 25% | 7.736937 | 1.000000e-09 | 0.998079 | 16.140808 | -11.0 | -3.125669 | -3.185303 | 6.016124 | 0.089999 | -0.532107 | -183.543793 | -0.519669 | -185.146241 | 0.0 |
| 50% | 18.351462 | 1.000000e-09 | 0.999633 | 36.912915 | -11.0 | 0.038935 | 0.003529 | 16.531034 | 0.092387 | -0.004279 | 1.527649 | 0.002749 | 0.106735 | 0.0 |
| 75% | 44.640179 | 1.000000e-09 | 0.999936 | 88.358651 | -11.0 | 3.224203 | 3.175715 | 43.416008 | 0.094875 | 0.518339 | 188.208423 | 0.525416 | 185.503679 | 0.0 |
| max | 4972.376090 | 1.000000e-09 | 1.000000 | 9731.697275 | -11.0 | 269.448714 | 82.801609 | 4972.880196 | 1.116313 | 16.176620 | 1570.451704 | 18.212887 | 1570.104422 | 0.0 |
plotSets = ['TransvPlane', 'TransvPsAngles', 'TransvPsMomenta', 'LongPsT']
plotDefs = bd.set_plot_defs_from_distrs([targetDistr], setNames=plotSets)
plotDefs[6]['lims2'] = [0., 100.]
_ = bd.plot_distr([targetDistr], plotDefs)
emitn = {}
emitGeom = {}
emitTraceSpace = {}
alphaTwiss = {}
betaTwiss = {}
gammaTwiss = {}
for planeName in ('x', 'y'):
emitn[planeName] = bd.compute_emittance(targetDistr, planeName, norm='normalized')
emitGeom[planeName] = bd.compute_emittance(targetDistr, planeName, norm='geometric')
emitTraceSpace[planeName] = bd.compute_emittance(targetDistr, planeName, norm='tracespace')
alphaTwiss[planeName], betaTwiss[planeName], gammaTwiss[planeName] = bd.compute_twiss(targetDistr, planeName)
print(
'emitn_{0:s} = {1:.1f} pi mm mrad, emitGeom_{0:s} = {2:.1f} pi mm mrad, emitTraceSpace_{0:s} = {3:.1f} pi mm mrad.'.format(
planeName, emitn[planeName], emitGeom[planeName], emitTraceSpace[planeName]))
print(
'alphaTwiss_{0:s} = {1:.3f}, betaTwiss_{0:s} = {2:.3f} mm, gammaTwiss_{0:s} = {3:.3f} 1/mm.'.format(
planeName, alphaTwiss[planeName], betaTwiss[planeName], gammaTwiss[planeName]))
Correcting offsets xAvg = -0.008 mm and pxAvg = 0.058 MeV/c. Correcting offsets xAvg = -0.008 mm and pxAvg = 0.058 MeV/c. Correcting offsets xAvg = -0.008 mm and xpAvg = 2.116 mrad. Correcting offsets xAvg = -0.008 mm and xpAvg = 2.116 mrad. Correcting offsets xAvg = -0.008 mm and xpAvg = 2.116 mrad. emitn_x = 14981.7 pi mm mrad, emitGeom_x = 153.8 pi mm mrad, emitTraceSpace_x = 484.8 pi mm mrad. alphaTwiss_x = -0.244, betaTwiss_x = 0.003 mm, gammaTwiss_x = 414.832 1/mm. Correcting offsets yAvg = 0.002 mm and pyAvg = 0.002 MeV/c. Correcting offsets yAvg = 0.002 mm and pyAvg = 0.002 MeV/c. Correcting offsets yAvg = 0.002 mm and ypAvg = 0.336 mrad. Correcting offsets yAvg = 0.002 mm and ypAvg = 0.336 mrad. Correcting offsets yAvg = 0.002 mm and ypAvg = 0.336 mrad. emitn_y = 14876.1 pi mm mrad, emitGeom_y = 152.7 pi mm mrad, emitTraceSpace_y = 484.3 pi mm mrad. alphaTwiss_y = -0.240, betaTwiss_y = 0.003 mm, gammaTwiss_y = 413.096 1/mm.
emitRef = emitTraceSpace['x']
gaussianPortions = {1: 0.6827, 2: 0.9545, 3: 0.9973, 4: 0.999937, 6: 0.99999998}
ellipseSpecs = {
'x': {'alphaTwiss': alphaTwiss['x'], 'betaTwiss': betaTwiss['x']},
'y': {'alphaTwiss': alphaTwiss['y'], 'betaTwiss': betaTwiss['y']}}
for planes in (['x'], ['x', 'y']):
print('\nFiltering on {:d} planes:'.format(len(planes)))
for Fa in (1., 2., 3., 4., 6.):
selEllipseSpecs = {k: v for k, v in ellipseSpecs.items() if k in planes}
distrWithinFaSigma, portionWithinFaSigma = bd.distr_within_ellipse(targetDistr, Fa**2.*emitRef, selEllipseSpecs)
print(
'Portion within {:.1f} sigma: {:.3f} (vs. {:.3f} for Gaussian).'.format(Fa, portionWithinFaSigma, gaussianPortions[Fa]))
Filtering on 1 planes: Portion within 1.0 sigma: 0.596 (vs. 0.683 for Gaussian). Portion within 2.0 sigma: 0.862 (vs. 0.955 for Gaussian). Portion within 3.0 sigma: 0.958 (vs. 0.997 for Gaussian). Portion within 4.0 sigma: 0.989 (vs. 1.000 for Gaussian). Portion within 6.0 sigma: 0.998 (vs. 1.000 for Gaussian). Filtering on 2 planes: Portion within 1.0 sigma: 0.448 (vs. 0.683 for Gaussian). Portion within 2.0 sigma: 0.793 (vs. 0.955 for Gaussian). Portion within 3.0 sigma: 0.932 (vs. 0.997 for Gaussian). Portion within 4.0 sigma: 0.980 (vs. 1.000 for Gaussian). Portion within 6.0 sigma: 0.996 (vs. 1.000 for Gaussian).
refEllipseSpecs = {
'x': {'alphaTwiss': alphaTwiss['x'], 'betaTwiss': betaTwiss['x']},
'y': {'alphaTwiss': alphaTwiss['y'], 'betaTwiss': betaTwiss['y']}
}
FaList = np.arange(6., 0, -0.2)
distrWithinRef = []
bunchPortions = []
emitnPortions = []
for Fa in FaList:
distrWithinRef.append(bd.distr_within_ellipse(targetDistr, Fa**2.*emitRef, refEllipseSpecs)[0])
bunchPortions.append(distrWithinRef[-1].shape[0]/targetDistr.shape[0])
emitnPortions.append(bd.compute_emittance(distrWithinRef[-1], 'x', verbose=False))
indsToAnnotate = [-5, -10, -15, -20, -25]
fig, ax = plt.subplots()
ax.plot(bunchPortions, emitnPortions, 'o-')
for ind in indsToAnnotate:
ax.annotate(np.round(FaList[ind], 1), (bunchPortions[ind]-0.05, emitnPortions[ind]))
ax.set_xlim([0, 1])
ax.set_ylim([0, 14e3])
ax.set_xlabel('Bunch portion')
ax.set_ylabel('Norm. transv. emittance [mm mrad]')
ax.grid()
# xy = np.column_stack([bunchPortions, emitnPortions])
# np.savetxt('emitnVsBunchPortion.dat', xy, header='{:24s} {:24s}'.format('bunchPortion', 'emitn [mm mrad]'))
fig, ax = plt.subplots()
ax.plot(FaList, bunchPortions, 'o-')
ax.set_xlim([0, np.max(FaList)])
ax.set_ylim([0, np.max(bunchPortions)])
ax.set_xlabel('Fa of cut with Fa^2*rmsEmit ellipse')
ax.set_ylabel('Bunch portion')
ax.grid()
stepsToPlot = np.arange(5*2, len(distrWithinRef), 2)
FaListToPlot = [FaList[stepInd] for stepInd in stepsToPlot]
plotDefs[1]['lims1'] = [-15., 15.]
plotDefs[2]['lims1'] = plotDefs[1]['lims1']
plotDefs[1]['displayTable'] = False
ax = bd.plot_distr(
[targetDistr, *[distrWithinRef[stepInd] for stepInd in stepsToPlot]], plotDefs,
legendLabels=[
'Full bunch',
*['Within {:.1f} sigma'.format(Fa) for Fa in FaListToPlot]])
for ind, Fa in enumerate(FaListToPlot):
bd.plot_ellipse(
ax[1][0,0], Fa**2.*emitRef, semiAxisOrder=1, alphaTwiss=alphaTwiss['x'], betaTwiss=betaTwiss['x'], color='k')
bd.plot_ellipse(
ax[2][0,0], Fa**2.*emitRef, semiAxisOrder=1, alphaTwiss=alphaTwiss['y'], betaTwiss=betaTwiss['y'], color='k')
EkinWindowEdges = np.arange(0., 500., 2.)
particlesInWindow = []
emitnWindows = []
for wInd in range(len(EkinWindowEdges)-1):
EkinWindow = {"Ekin": EkinWindowEdges[wInd:wInd+2]}
particlesInWindow.append(bd.filter_distr(targetDistr, EkinWindow))
emitnWindows.append(bd.compute_emittance(particlesInWindow[wInd], 'x', verbose=False))
fig, ax = plt.subplots()
ax.plot(EkinWindowEdges[:-1], emitnWindows, 'o-')
ax.set_ylim([0, 15e3])
ax.set_xlabel('Ekin window bottom edge [MeV]')
ax.set_ylabel('Norm. emittance [pimmmrad]')
ax.grid()
windowIndsToPlot = np.arange(0, len(particlesInWindow), 1)
windowsToPlot = [EkinWindowEdges[wInd:wInd+2] for wInd in windowIndsToPlot]
ax = bd.plot_distr(
[targetDistr, *[particlesInWindow[wInd] for wInd in windowIndsToPlot]], plotDefs,
legendLabels=[
'Full bunch',
*['{:.1f} MeV < Ekin < {:.1f} MeV'.format(*EkinWindow) for EkinWindow in windowsToPlot]])
/home/tia/Repos/GIT_PSIPositronProduction/BdPyTools/BeamDynamics.py:1123: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`). Consider using `matplotlib.pyplot.close()`. fig, ax = plt.subplots(2, 2, figsize=(figWidth, figHeight))
DISTR_REL_PATH = {}
FILTER_SPECS = {}
REF_PARTICLE_PZ = {}
DISTR_REL_PATH['wp1'] = '../../Data/Astra/WP1atDR.dat'
FILTER_SPECS['wp1'] = {
'z': [14480., 14530.],
'pz': [0., 1620.]}
REF_PARTICLE_PZ['wp1'] = 1565. # [MeV/c]
DISTR_REL_PATH['wp2'] = '../../Data/Astra/WP2atDR.dat'
FILTER_SPECS['wp2'] = {
'z': [14470., 14520.],
'pz': [0., 1620.]}
REF_PARTICLE_PZ['wp2'] = 1565. # [MeV/c]
drDistr = {}
for wpName, relPath in DISTR_REL_PATH.items():
drDistr[wpName] = pd.read_csv(relPath, names=bd.FILE_TYPES_SPECS['astra']['columnOrder'], header=None)
# Minimal substitution of drDistr[wpName] = bd.extend_standard_df(drDistr[wpName], False)
drDistr[wpName]['xp'] = 0.
drDistr[wpName]['yp'] = 0.
drDistr[wpName]['t'] = 0.
drDistr[wpName]['Ekin'] = 0.
drDistr[wpName]['pdgId'] = -11
print('Working point {:s}:'.format(wpName))
print(drDistr[wpName].describe())
Working point wp1:
x y z px py pz \
count 30451.0 30451.0 30451.000000 30451.0 30451.0 30451.000000
mean 0.0 0.0 14451.635141 0.0 0.0 1308.134190
std 0.0 0.0 605.428368 0.0 0.0 422.126061
min 0.0 0.0 -9991.133519 0.0 0.0 0.567927
25% 0.0 0.0 14494.416345 0.0 0.0 1183.249370
50% 0.0 0.0 14506.803113 0.0 0.0 1533.458841
75% 0.0 0.0 14511.485757 0.0 0.0 1575.519662
max 0.0 0.0 64988.153966 0.0 0.0 4883.870272
clock macroCharge particleIndex statusFlag xp yp \
count 30451.0 30451.0 30451.0 30451.0 30451.0 30451.0
mean 0.0 0.0 0.0 0.0 0.0 0.0
std 0.0 0.0 0.0 0.0 0.0 0.0
min 0.0 0.0 0.0 0.0 0.0 0.0
25% 0.0 0.0 0.0 0.0 0.0 0.0
50% 0.0 0.0 0.0 0.0 0.0 0.0
75% 0.0 0.0 0.0 0.0 0.0 0.0
max 0.0 0.0 0.0 0.0 0.0 0.0
t Ekin pdgId
count 30451.0 30451.0 30451.0
mean 0.0 0.0 -11.0
std 0.0 0.0 0.0
min 0.0 0.0 -11.0
25% 0.0 0.0 -11.0
50% 0.0 0.0 -11.0
75% 0.0 0.0 -11.0
max 0.0 0.0 -11.0
Working point wp2:
x y z px py pz \
count 29279.0 29279.0 29279.000000 29279.0 29279.0 29279.000000
mean 0.0 0.0 14459.170166 0.0 0.0 1421.584427
std 0.0 0.0 366.967385 0.0 0.0 335.120112
min 0.0 0.0 -11699.220033 0.0 0.0 1.290207
25% 0.0 0.0 14492.560951 0.0 0.0 1501.087186
50% 0.0 0.0 14496.303136 0.0 0.0 1557.978180
75% 0.0 0.0 14497.844773 0.0 0.0 1574.104075
max 0.0 0.0 21798.301070 0.0 0.0 3536.281605
clock macroCharge particleIndex statusFlag xp yp \
count 29279.0 29279.0 29279.0 29279.0 29279.0 29279.0
mean 0.0 0.0 0.0 0.0 0.0 0.0
std 0.0 0.0 0.0 0.0 0.0 0.0
min 0.0 0.0 0.0 0.0 0.0 0.0
25% 0.0 0.0 0.0 0.0 0.0 0.0
50% 0.0 0.0 0.0 0.0 0.0 0.0
75% 0.0 0.0 0.0 0.0 0.0 0.0
max 0.0 0.0 0.0 0.0 0.0 0.0
t Ekin pdgId
count 29279.0 29279.0 29279.0
mean 0.0 0.0 -11.0
std 0.0 0.0 0.0
min 0.0 0.0 -11.0
25% 0.0 0.0 -11.0
50% 0.0 0.0 -11.0
75% 0.0 0.0 -11.0
max 0.0 0.0 -11.0
for wpName, filterSpecs in zip(drDistr.keys(), FILTER_SPECS.values()):
drDistr[wpName] = bd.filter_distr(drDistr[wpName], filterSpecs)
print('Working point {:s}:'.format(wpName))
print(drDistr[wpName].describe())
Working point wp1:
x y z px py pz \
count 25337.0 25337.0 25337.000000 25337.0 25337.0 25337.000000
mean 0.0 0.0 14506.027271 0.0 0.0 1414.660221
std 0.0 0.0 7.661450 0.0 0.0 321.748731
min 0.0 0.0 14484.702394 0.0 0.0 5.577844
25% 0.0 0.0 14501.894891 0.0 0.0 1458.721448
50% 0.0 0.0 14508.587808 0.0 0.0 1549.095748
75% 0.0 0.0 14511.954643 0.0 0.0 1579.218798
max 0.0 0.0 14527.074149 0.0 0.0 1619.837670
clock macroCharge particleIndex statusFlag xp yp \
count 25337.0 25337.0 25337.0 25337.0 25337.0 25337.0
mean 0.0 0.0 0.0 0.0 0.0 0.0
std 0.0 0.0 0.0 0.0 0.0 0.0
min 0.0 0.0 0.0 0.0 0.0 0.0
25% 0.0 0.0 0.0 0.0 0.0 0.0
50% 0.0 0.0 0.0 0.0 0.0 0.0
75% 0.0 0.0 0.0 0.0 0.0 0.0
max 0.0 0.0 0.0 0.0 0.0 0.0
t Ekin pdgId
count 25337.0 25337.0 25337.0
mean 0.0 0.0 -11.0
std 0.0 0.0 0.0
min 0.0 0.0 -11.0
25% 0.0 0.0 -11.0
50% 0.0 0.0 -11.0
75% 0.0 0.0 -11.0
max 0.0 0.0 -11.0
Working point wp2:
x y z px py pz \
count 26491.0 26491.0 26491.000000 26491.0 26491.0 26491.000000
mean 0.0 0.0 14494.991173 0.0 0.0 1489.906379
std 0.0 0.0 5.332383 0.0 0.0 232.064203
min 0.0 0.0 14472.672312 0.0 0.0 8.290856
25% 0.0 0.0 14494.377212 0.0 0.0 1528.886562
50% 0.0 0.0 14496.620400 0.0 0.0 1561.333039
75% 0.0 0.0 14497.979419 0.0 0.0 1575.815748
max 0.0 0.0 14519.997744 0.0 0.0 1598.811507
clock macroCharge particleIndex statusFlag xp yp \
count 26491.0 26491.0 26491.0 26491.0 26491.0 26491.0
mean 0.0 0.0 0.0 0.0 0.0 0.0
std 0.0 0.0 0.0 0.0 0.0 0.0
min 0.0 0.0 0.0 0.0 0.0 0.0
25% 0.0 0.0 0.0 0.0 0.0 0.0
50% 0.0 0.0 0.0 0.0 0.0 0.0
75% 0.0 0.0 0.0 0.0 0.0 0.0
max 0.0 0.0 0.0 0.0 0.0 0.0
t Ekin pdgId
count 26491.0 26491.0 26491.0
mean 0.0 0.0 -11.0
std 0.0 0.0 0.0
min 0.0 0.0 -11.0
25% 0.0 0.0 -11.0
50% 0.0 0.0 -11.0
75% 0.0 0.0 -11.0
max 0.0 0.0 -11.0
plotSets = ['z-pz']
plotDefs = bd.set_plot_defs_from_distrs(drDistr.values(), setNames=plotSets)
plotDefs[0]['lims1'] = [14488., 14516.]
plotDefs[0]['lims2'] = [1420., 1620.]
for distr in drDistr.values():
_ = bd.plot_distr([distr], plotDefs)
def particles_in_window(distr, refParticlePz, pzWindowHalfAmplitudes, pzWindowTopEdgeScanRange, pzWindowAmplitudeStep):
particlesInWindow = []
for pzWindowAmplitude in pzWindowHalfAmplitudes:
optYield = 0
PZ_WINDOW_TOP_STEP = 1. # [MeV/c]
for pzWindowTop in np.arange(distr['pz'].max(), distr['pz'].max()-pzWindowTopEdgeScanRange*(1+PZ_WINDOW_TOP_STEP), -PZ_WINDOW_TOP_STEP):
pzWindow = {
"pz": [pzWindowTop-2.*pzWindowAmplitude*refParticlePz, pzWindowTop]
}
tmpParticlesInWindow = bd.filter_distr(distr, pzWindow)
if tmpParticlesInWindow.shape[0] > optYield:
optYield = tmpParticlesInWindow.shape[0]
optParticlesInWindow = tmpParticlesInWindow
particlesInWindow.append(optParticlesInWindow)
return particlesInWindow
REF_WINDOW_HALF_AMPLITUDE = 0.038
MAX_WINDOW_HALF_AMPLITUDE = 0.06
PZ_WINDOW_AMPLITUDE_STEP = 0.001
PZ_WINDOW_TOP_EDGE_SCAN_RANGE = 200. # [MeV/c]
pzWindowHalfAmplitudes = np.flip(np.arange(PZ_WINDOW_AMPLITUDE_STEP, MAX_WINDOW_HALF_AMPLITUDE+PZ_WINDOW_AMPLITUDE_STEP, PZ_WINDOW_AMPLITUDE_STEP))
particlesInWindow = {}
totParticlesInWindow = {}
for wpName in drDistr.keys():
particlesInWindow[wpName] = particles_in_window(drDistr[wpName], REF_PARTICLE_PZ[wpName], pzWindowHalfAmplitudes, PZ_WINDOW_TOP_EDGE_SCAN_RANGE, PZ_WINDOW_AMPLITUDE_STEP)
totParticlesInWindow[wpName] = [partInWin.shape[0] / drDistr[wpName].shape[0] for partInWin in particlesInWindow[wpName]]
totParticlesInRefWindow = np.interp(REF_WINDOW_HALF_AMPLITUDE, np.flip(pzWindowHalfAmplitudes), np.flip(totParticlesInWindow[wpName]))
totParticlesInWindow[wpName] /= totParticlesInRefWindow
fig, ax = plt.subplots()
for wpName in drDistr.keys():
ax.plot(pzWindowHalfAmplitudes*100., totParticlesInWindow[wpName], '-')
ax.plot(REF_WINDOW_HALF_AMPLITUDE*100., 1., 'ro')
ax.set_xlim([0., 5.])
ax.set_ylim([0., 1.2])
ax.set_xlabel('pz window amplitude [%]')
ax.set_ylabel('Particles in window')
ax.legend([wpName.upper() for wpName in drDistr.keys() for _ in range(2)])
ax.grid()
# xy1y2 = np.column_stack([pzWindowHalfAmplitudes*100., *totParticlesInWindow.values()])
# np.savetxt('drAcceptanceVsPzWindowAmplitude.dat', xy1y2, header='{:24s} {:24s} {:24s}'.format('pzWindowHalfAmpl [%]', 'DRAcceptanceWP1', 'DRAcceptanceWP2'))
stepsToPlot = np.arange(0, len(particlesInWindow['wp1']), 10)
pzWindowHalfAmplitudesToPlot = [pzWindowHalfAmplitudes[stepInd] for stepInd in stepsToPlot]
for wpName in drDistr:
ax = bd.plot_distr(
[drDistr[wpName], *[particlesInWindow[wpName][stepInd] for stepInd in stepsToPlot]], plotDefs,
legendLabels=[
'Full bunch',
*['Within +- {:.1f} %'.format(pzWindowHalfAmplitude*100.) for pzWindowHalfAmplitude in pzWindowHalfAmplitudesToPlot]])
